library(MendelianRandomization)
Registered S3 method overwritten by 'dplyr':
method from
print.rowwise_df
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
data_dir = "../../dat"
seq_predictions_all = read.csv(file.path(data_dir, "means_and_uncertainties.csv"))
seq_predictions = subset(seq_predictions_all, seq_num == 1)
str(seq_predictions)
'data.frame': 3001 obs. of 6 variables:
$ seq_num : int 1 1 1 1 1 1 1 1 1 1 ...
$ mut : Factor w/ 3000 levels "0A","0C","0G",..: 2729 1 2 3 342 343 344 674 675 676 ...
$ X_pred_mean: num 2.33e-03 3.50e-05 -1.89e-05 -1.57e-05 -1.33e-05 ...
$ X_pred_var : num 0.0586 0.06 0.06 0.06 0.06 ...
$ Y_pred_mean: num 1.71e-04 8.74e-06 -2.57e-06 -1.09e-06 -2.20e-06 ...
$ Y_pred_var : num 0.00651 0.00777 0.00777 0.00777 0.00777 ...
bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])
MRInputObject <- mr_input(bx = bx,
bxse = bxse,
by = by,
byse = byse)
cis = c()
for (seq in (1:50)) {
seq_predictions = subset(seq_predictions_all, seq_num == seq)
bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])
MRInputObject <- mr_input(bx = bx,
bxse = bxse,
by = by,
byse = byse)
EggerObject <- mr_egger(MRInputObject, robust = FALSE, alpha = .05)
# cis = append(cis, c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
print(c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
}
[1] 0.0590192 0.1159848
[1] 0.3233909 0.5597846
[1] 0.03235503 0.05877458
[1] 0.1188956 0.2318077
[1] 0.1645073 0.2648561
[1] 0.2289623 0.4052021
[1] 0.3443012 0.6641439
[1] 0.1756041 0.3033061
[1] 0.1583441 0.2399981
[1] 0.1889218 0.3005885
[1] 0.1827082 0.3066469
[1] 0.05220543 0.10226116
[1] 0.91425 1.57727
[1] 0.2863903 0.4708377
[1] 0.1051252 0.1822457
[1] 0.03974643 0.07074950
[1] 0.07505532 0.11982899
[1] 0.3706091 0.5714237
[1] 0.4507038 0.7233038
[1] 0.2518358 0.4329524
[1] 0.02264022 0.03560069
[1] 0.2323919 0.3507465
[1] 0.07638818 0.12772401
[1] 0.3171318 0.4806720
[1] 0.2739240 0.4173557
[1] 0.03117739 0.06000044
[1] 0.04368776 0.12557530
[1] 0.07505245 0.14740756
[1] 0.04114612 0.11159205
[1] 0.1269643 0.2866314
[1] 0.2601743 0.6519868
[1] 0.06622755 0.11109060
[1] 1.109005 1.907397
[1] 0.04254034 0.09048482
[1] 0.1153641 0.1847817
[1] 0.4348032 0.9299903
[1] 0.1471299 0.3256183
[1] 0.08097467 0.16127508
[1] 0.2559104 0.4957051
[1] 0.4693041 0.8698878
[1] 0.02719744 0.05512086
[1] 0.1165405 0.2815027
[1] 0.1989373 0.4144706
[1] 0.06628877 0.14424319
[1] 0.1635376 0.4410249
[1] 0.1459457 0.3396873
[1] 0.03371326 0.09627066
[1] 0.1937338 0.4058515
[1] 0.9633243 1.6609768
[1] 0.1370262 0.3841660
# dim(cis) <- c(2, 50)
# print(cis)
mr_plot(MRInputObject)
install.packages("devtools")
Installing package into ‘/home/stephenmalina/R/x86_64-pc-linux-gnu-library/3.6’
(as ‘lib’ is unspecified)
trying URL 'https://cloud.r-project.org/src/contrib/devtools_2.2.1.tar.gz'
Content type 'application/x-gzip' length 372273 bytes (363 KB)
==================================================
downloaded 363 KB
* installing *source* package ‘devtools’ ...
** package ‘devtools’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
*** copying figures
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (devtools)
The downloaded source packages are in
‘/tmp/RtmpfWDaCQ/downloaded_packages’
library(devtools)
Loading required package: usethis
install_github("WSpiller/RadialMR")
Skipping install of 'RadialMR' from a github remote, the SHA1 (18a7a6e9) has not changed since last install.
Use `force = TRUE` to force installation
plot(by/bx, bx/byse, xlim=c(min((by-2*byse)/bx), max((bx+2*byse)/bx)), ylim=c(0, max(bx/byse)))
for (j in 1:length(bx)) {
lines(c((by[j]-1.96*byse[j])/bx[j], (by[j]+1.96*byse[j])/bx[j]),
c(bx[j]/byse[j], bx[j]/byse[j]))
}
abline(h=0, lwd=1); abline(v=0, lwd=1)

library(RadialMR)
r_input <- format_radial(bx, by, bxse, byse)
Missing SNP IDs; Generating placeholders
funnel_radial(egger_radial(r_input, alpha = .05))
Want to understand how all the pieces fit together? See the R for Data Science book: http://r4ds.had.co.nz/
Weights not specified: Adopting modified second-order weights
Radial MR-Egger
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01214285 0.014905603 -0.8146503 0.4153372
Wj 0.44241552 0.005847025 75.6650678 0.0000000
Residual standard error: 0.728 on 2999 degrees of freedom
F-statistic: 5725.2 on 1 and 2999 DF, p-value: <2e-16
Q-Statistic for heterogeneity: 1588.522 on 2999 DF , p-value: 1
Outliers detected

LS0tCnRpdGxlOiAiTWVuZGVsaWFuIFJhbmRvbWl6YXRpb24gYW5hbHlzaXMgb2YgRGVlcFNFQSBnZW5lcmF0ZWQgZGF0YSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KbGlicmFyeShNZW5kZWxpYW5SYW5kb21pemF0aW9uKQpgYGAKCmBgYHtyfQpkYXRhX2RpciA9ICIuLi8uLi9kYXQiCnNlcV9wcmVkaWN0aW9uc19hbGwgPSByZWFkLmNzdihmaWxlLnBhdGgoZGF0YV9kaXIsICJtZWFuc19hbmRfdW5jZXJ0YWludGllcy5jc3YiKSkKc2VxX3ByZWRpY3Rpb25zID0gc3Vic2V0KHNlcV9wcmVkaWN0aW9uc19hbGwsIHNlcV9udW0gPT0gMSkKYGBgCgpgYGB7cn0Kc3RyKHNlcV9wcmVkaWN0aW9ucykKYGBgCmBgYHtyfQpieCA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJYX3ByZWRfbWVhbiJdKQpieHNlIDwtIHVubGlzdChzZXFfcHJlZGljdGlvbnNbIlhfcHJlZF92YXIiXSkKYnkgPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWV9wcmVkX21lYW4iXSkKYnlzZSA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJZX3ByZWRfdmFyIl0pCmBgYAoKYGBge3J9Ck1SSW5wdXRPYmplY3QgPC0gbXJfaW5wdXQoYnggPSBieCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgYnhzZSA9IGJ4c2UsIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ5ID0gYnksIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ5c2UgPSBieXNlKQpgYGAKCgpgYGB7cn0KY2lzID0gYygpCmZvciAoc2VxIGluICgxOjUwKSkgewogIHNlcV9wcmVkaWN0aW9ucyA9IHN1YnNldChzZXFfcHJlZGljdGlvbnNfYWxsLCBzZXFfbnVtID09IHNlcSkKICAKICBieCA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJYX3ByZWRfbWVhbiJdKQogIGJ4c2UgPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWF9wcmVkX3ZhciJdKQogIGJ5IDwtIHVubGlzdChzZXFfcHJlZGljdGlvbnNbIllfcHJlZF9tZWFuIl0pCiAgYnlzZSA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJZX3ByZWRfdmFyIl0pCiAgTVJJbnB1dE9iamVjdCA8LSBtcl9pbnB1dChieCA9IGJ4LCAKICAgICAgICAgICAgICAgICAgICAgICAgICBieHNlID0gYnhzZSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgYnkgPSBieSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgYnlzZSA9IGJ5c2UpCiAgRWdnZXJPYmplY3QgPC0gbXJfZWdnZXIoTVJJbnB1dE9iamVjdCwgcm9idXN0ID0gRkFMU0UsIGFscGhhID0gLjA1KQogICMgY2lzID0gYXBwZW5kKGNpcywgYyhFZ2dlck9iamVjdCRDSUxvd2VyLkVzdCwgRWdnZXJPYmplY3QkQ0lVcHBlci5Fc3QpKQogIHByaW50KGMoRWdnZXJPYmplY3QkQ0lMb3dlci5Fc3QsIEVnZ2VyT2JqZWN0JENJVXBwZXIuRXN0KSkKfQojIGRpbShjaXMpIDwtIGMoMiwgNTApCiMgcHJpbnQoY2lzKQpgYGAKCgpgYGB7cn0KbXJfcGxvdChNUklucHV0T2JqZWN0KQpgYGAKYGBge3J9Cmluc3RhbGwucGFja2FnZXMoImRldnRvb2xzIikKYGBgCgpgYGB7cn0KbGlicmFyeShkZXZ0b29scykKaW5zdGFsbF9naXRodWIoIldTcGlsbGVyL1JhZGlhbE1SIikKYGBgCmBgYHtyfQpwbG90KGJ5L2J4LCBieC9ieXNlLCB4bGltPWMobWluKChieS0yKmJ5c2UpL2J4KSwgbWF4KChieCsyKmJ5c2UpL2J4KSksIHlsaW09YygwLCBtYXgoYngvYnlzZSkpKQpmb3IgKGogaW4gMTpsZW5ndGgoYngpKSB7CmxpbmVzKGMoKGJ5W2pdLTEuOTYqYnlzZVtqXSkvYnhbal0sIChieVtqXSsxLjk2KmJ5c2Vbal0pL2J4W2pdKSwKYyhieFtqXS9ieXNlW2pdLCBieFtqXS9ieXNlW2pdKSkKfQphYmxpbmUoaD0wLCBsd2Q9MSk7IGFibGluZSh2PTAsIGx3ZD0xKQpgYGAKYGBge3J9CmxpYnJhcnkoUmFkaWFsTVIpCnJfaW5wdXQgPC0gZm9ybWF0X3JhZGlhbChieCwgYnksIGJ4c2UsIGJ5c2UpCmZ1bm5lbF9yYWRpYWwoZWdnZXJfcmFkaWFsKHJfaW5wdXQsIGFscGhhID0gLjA1KSkKYGBgCgo=